A scenario modelling analysis to anticipate the impact of COVID-19 vaccination in adolescents and children on disease outcomes in the Netherlands, summer 2021

Background Since the roll-out of COVID-19 vaccines in late 2020 and throughout 2021, European governments have relied on mathematical modelling to inform policy decisions about COVID-19 vaccination. Aim We present a scenario-based modelling analysis in the Netherlands during summer 2021, to inform whether to extend vaccination to adolescents (12–17-year-olds) and children (5–11-year-olds). Methods We developed a deterministic, age-structured susceptible-exposed-infectious-recovered (SEIR) model and compared modelled incidences of infections, hospital and intensive care admissions, and deaths per 100,000 people across vaccination scenarios, before the emergence of the Omicron variant. Results Our model projections showed that, on average, upon the release of all non-pharmaceutical control measures on 1 November 2021, a large COVID-19 wave may occur in winter 2021/22, followed by a smaller, second wave in spring 2022, regardless of the vaccination scenario. The model projected reductions in infections/severe disease outcomes when vaccination was extended to adolescents and further reductions when vaccination was extended to all people over 5 years-old. When examining projected disease outcomes by age group, individuals benefitting most from extending vaccination were adolescents and children themselves. We also observed reductions in disease outcomes in older age groups, particularly of parent age (30–49 years), when children and adolescents were vaccinated, suggesting some prevention of onward transmission from younger to older age groups. Conclusions While our scenarios could not anticipate the emergence/consequences of SARS-CoV-2 Omicron variant, we illustrate how our approach can assist decision making. This could be useful when considering to provide booster doses or intervening against future infection waves.


Methods
We developed a deterministic age-structured compartmental susceptible-exposedinfectious-recovered (SEIR) model extended to include states for severe disease outcomes and vaccination status. The population is partitioned into 10-year age groups (0-9, 10-19, …, 70-79, 80+). Within each age group we further partition the population into those who are unvaccinated, and separately for those who are vaccinated with 1 to 5 doses, and then finally into disease states: susceptible (S), infected but not yet infectious (E), infectious (I), hospitalized (H), in intensive care (IC), return to the hospital ward after intensive care (HIC), recovered (R), and dead (D) ( Figure S1).
When a person is vaccinated, they first enter a hold state where they are vaccinated, but not yet (fully) protected. After a delay period, they enter the vaccinated and protected state for the dose they have received. We assume only susceptible individuals are vaccinated. We use the model to determine the number of daily cases, hospital admissions, and IC admissions under different vaccination scenarios. The mathematical equations for determining each outcome and the differential model equations are shown below. In the model equations, bold capital letters refer to matrices, bold lower case letters indicate vectors, and plain text lower case letters indicate scalars. Since model compartments are denoted with capital letters, but are vectors due to the agestratification of the model, we denote them with a half arrow above the compartment symbol, e.g., ⃑ . Parameter definitions and values are shown in Table S3 and Table S4.

Outcome Equations
The model is designed to incorporate a single vaccine product with a 2-dose regimen that 1) reduces susceptibility to infection, 2) reduces risk of hospitalization if a vaccinated individual is infected, and 3) reduces risk of infecting others (transmission) if a vaccinated person is infected. However, since there are more than one vaccine products currently licensed for use against COVID-19 in The Netherlands (the vaccines made by Pfizer/BioNTech, Moderna, AstraZeneca, and Janssen), we incorporate different vaccine products by taking the daily weighted average of the number of people with each vaccine product (and dose), the corresponding delay to protection of each vaccine product, and the vaccine effectiveness against each outcome (Table S1). Janssen is incorporated by using zero for the number of second doses at all time points. The weighted average of the VE can be expressed as: where ( ) is the number of vaccines given as dose d (d = 1, 2) from vaccine product i (i = 1, 2, 3, 4) at time t and is the VE against outcome o (o = infection, hospitalization, transmission), for dose d, from vaccine product i.
When VE is assumed to wane, we include the amount of waning as a logistic function parameterized so that after 6 months vaccine-induced protection is reduced by 50% and reduced by 100% after 1 year. The amount of waning is calculated at a given time since vaccination as where k is the logistic growth rate (here, k = 0.03), t is the time since vaccination in days, and t0 is the time point (in days) where 50% reduction occurs (here, t0=180). Then at each timepoint, VE with waning ( ( )) is calculated as Rate of vaccination for each day, vaccine product, and dose for each age group is a model input. Increase of vaccination coverage over time was based on weekly data of allocated vaccines (up to 16 June 2021) or projected available vaccines (after 16 June 2021) by vaccine type, dose number, and target group (split up between healthcare workers, residents in institutions, chronically ill, and, if not part of one of the aforementioned groups, by age in 10-years age-bands) [1]. Projected final vaccination coverages were assumed at 90% for ≥70-year-olds and residents in institutions, 85% for 60-to 69-years-olds, health care workers and high-risk individuals, and 75% for others below 60 years of age. It was assumed that available vaccines were administered immediately until a target group reached the final coverage.
To account for the seasonal pattern of SARS-CoV-2 whereby, transmission is lower in summer and higher in winter, we define the transmission rate at time t, ( ), as a sinusoidal function of seasonality [2]: 0 is a baseline (non-seasonal) transmission rate, 1 is the amplitude of seasonal forcing, and t is the day of the year. We assume 1 = 0.14. The baseline transmission rate 0 and initial conditions for forward simulations are estimated by fitting the model to daily cases from the national notification database Osiris from 01 January 2020 to 22 June 2021, when vaccination in 12-17 year olds began in The Netherlands ( Figure S2). The model is fitted to data piecewise to correspond with the correct contact patterns associated with different non-pharmaceutical interventions within each time window (Table S2). We directly estimate effective reproduction number (Rt) within each time window using maximum likelihood estimation. We assume daily cases follow a negative binomial distribution with mean µ and overdispersion parameter . Estimates of Rt are converted to transmission rate by: 0 = * , where γ is the inverse of the infectious period and ρ is the dominant eigen value of the product of the diagonal matrix of the number of susceptibles and the transmission matrix. The transmission matrix is determined by multiplying the rows and columns of the contact matrix by estimates of the relative susceptibility and infectiousness of each age group compared to the 0-9 years age group (Table S4).
The amplitude of seasonal forcing 1 is estimated by linear regression: The response variable is the logarithm of the effective reproduction number [3] divided by the fraction susceptible s(t). The fraction susceptible is estimated from the cumulative age-specific hospitalizations at time t divided by age-specific hospitalization rates, which were in turn estimated from seroprevalences after the first wave in The Netherlands [3]. The explanatory variables for seasonal changes are absolute humidity (XAH)and temperature in degrees Celsius (Xtemp) [4]. Additional explanatory variables that may affect were also included: percent change in mobility from "workplace" and "retail and recreation" (Xmobility) [5], increased transmissibility due to change in circulating virus variants (notably alpha and delta) (Xvariant) [6], day of the week (Xday), and intervention period (Xintervention). Intervention period is a variable indicating 2-week to 2month time periods during the epidemic without change in control measures. New periods started 2020-03-13, 2020-05-06, 2020-06-01, 2020-07-01, 2020-09-01, 2020-09-29, 2020-10-26, 2020-12-21, 2021-01-23, 2021-02-08, 2021-03-01, 2021-04-26, 2021-05-17, 2021-06-05. The seasonality curve was extracted from the intercept and seasonal variable coefficients only, from which the amplitude was determined by fitting a sinusoidal curve. a A delay of protection is assumed to be 14 days after the first dose and 7 days after the second dose for the Pfizer/BioNTech, Moderna, and AstraZeneca vaccines and 28 days after the first dose of the Janssen vaccine. b VE against hospitalization is incorporated as a multiplier on the probability of being hospitalized after infection as (1 -VEhospitalization)/(1-VEinfection) to account for the inclusion of people who are never infected (and thus never hospitalized) included in the estimation of VE against hospitalization.   Table S1 With multiple vaccines, the weighted average is used η 1 -vaccine efficacy of first dose See Table S1 With multiple vaccines, the weighted average is used 1 -vaccine efficacy of second dose See Table S1 With multiple vaccines, the weighted average is used